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Abstract The recently developed Flexible Local Approximation MEthod 
(FLAME) produces accurate difference schemes by replacing the usual Tay- 
lor expansion with Trefftz functions - local solutions of the underlying dif- 
ferential equation. This paper advances and casts in a general form a sig- 
nificant modification of FLAME proposed recently by Pinheiro & Webb: 
a least-squares fit instead of the exact match of the approximate solution 
at the stencil nodes. As a consequence of that, FLAME schemes can now 
be generated on irregular stencils with the number of nodes substantially 
greater than the number of approximating functions. The accuracy of the 
method is preserved but its robustness is improved. For demonstration, the 
paper presents a number of numerical examples in 2D and 3D: electrostatic 
(magnetostatic) particle interactions, scattering of electromagnetic (acous- 
tic) waves, and wave propagation in a photonic crystal. The examples ex- 
plore the role of the grid and stencil size, of the number of approximating 
functions, and of the irregularity of the stencils. 

Key words Flexible local approximation, finite difference schemes, Tr- 
efftz functions, electrostatics, wave propagation, wave scattering, electro- 
magnetic waves, multiparticle problems, long-range interactions, irregular 
stencils, meshless methods, least squares approximation. 



1 Introduction 



Traditional finite difference analysis relies primarily on Taylor expansions. 
These are quite general but are accurate only if the underlying solution has 
the level of smoothness commensurate with the order of the expansion. In 
particular, the Taylor approximation breaks down at material interfaces due 
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to jumps in the solution and/or its derivatives. This leads to the well known 
"staircase effect" in difference schemes (see a discussion of the algebraic na- 
ture of this numerical artifact in [1,2]). Other common situations where the 
Taylor series, and hence the corresponding classical schemes, are inaccurate 
include boundary layers and sharp peaks or singularities in the vicinity of 
sources, edges and corners. 

The Flexible Local Approximation MEthod (FLAME) [1,2] replaces the 
Taylor polynomials with much more accurate "Trefftz" approximations by 
local solutions of the underlying differential equation. As a result, the ap- 
proximation accuracy and consequently the consistency error of the scheme 
can be improved dramatically. 

Previously, FLAME relied on the point-matching of the nodal values of 
the local approximation. To fix the key idea, consider a nine-point (3 x 3) 
scheme for the Laplace equation in 2D. The solution is approximated lo- 
cally as a linear combination of eight Trefftz functions, harmonic polyno- 
mials being the most natural choice.^ With nine nodes and only eight free 
parameters, the nodal values must be linearly dependent. It is the linear re- 
lationship between them that constitutes the FLAME scheme. Details and 
a large variety of examples can be found in [1,2,4,5,6,7]; the applications 
include electro- and magnetistatics, wave propagation and scattering, the 
Poisson-Boltzmann equation in macromolecular and colloidal simulation. 

One limitation of the point matching procedure is that it links the num- 
ber of approximating functions and stencil nodes, as evident from the above 
example of eight basis functions for the nine-point stencil. This connection 
between functions and nodes has impeded further progress of the method 
in two different ways. 

First, the most natural choices of the stencil and the basis set have not 
always been feasible. For instance, consider the standard 3x3 stencil in 
2D problems involving cylindrical particles. A natural choice of the basis 
in such cases is cylindrical harmonics [2,5,7] that include the exponential 
factors cxp(im0), m = 0, ±1, ±2, . . ., where (j) is the polar angle. Since these 
harmonics come in pairs ±m for all m except m = 0, the "natural" number 
of functions in this basis is odd. At the same time, to obtain a FLAME 
scheme on the nine-point stencil, one must choose eight basis functions, 
which breaks the natural symmetry: only one index m from the pair m = ±4 
can be taken. This is more of a nuisance than an actual detriment because 
all harmonics up to order m = 3 are still included. 

A more serious obstacle is that large stencils (in terms of the number 
of nodes) require in the established version of FLAME a commensurately 
large number of basis functions. Unfortunately, expanded basis sets tend 
to be poorly conditioned. Yet large stencils are highly desirable in some 
circumstances, in particular for irregular distributions of nodes where small 



^ There are of course other possibilities: e.g. functions r™ exp(im(;/)) (m = 
0, ±1, ±2, . . .) in a polar coordinate system (r, 0) or Green's functions of the form 
|r — Tare I ^'^j where the source rare is located away from the given stencil [3], etc. 
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stencils can be strongly distorted and unreliable [8]; including more nodes 
in the stencil tends to increase robustness. 

The disparity between the desirable number of basis functions and stencil 
nodes became particularly apparent in three-dimensional electromagnetic 
vector problems. If the number of basis functions is commensurate with 
the number of the nodal degrees of freedom (number of stencil nodes times 
three Cartesian components), the basis sets become ill-conditioned. Pinheiro 
& Webb recently overcame this difficulty by using a relatively small number 
of approximating functions and applying a least-squares match of the nodal 
values [9]. It is this idea that is further advanced in the present paper. 

The use of least squares matching allows one to relax or even sever the 
connection between the basis functions and the stencil. This facilitates the 
construction of FLAME on non-canonical irregular stencils. So far FLAME 
has mostly been used on regular Cartesian grids - not as a requirement but 
as a practical matter, to avoid dealing with distorted / skewed stencils. An 
exception to this practice was adaptive FLAME in [8] , where complications 
due to irregular distributions of nodes in adaptive refinement stencils had 
to be overcome. 

Irregular and adaptive node arrangements are in general highly desirable 
because the nodes can be concentrated in areas where they are needed the 
most, in contrast with a regular grid that can only be refined globally. 
The least squares FLAME described in this paper works on non-canonical 
stencils more reliably than the previous versions of FLAME. 

Irregular stencils are clearly a feature that the new approach, least- 
squares FLAME, shares with meshless methods (see reviews [10,11]). An- 
other feature that is shared - albeit superficially (see below) - is the least 
squares matching. I shall, however, refrain from referring to least squares 
FLAME as a meshless method, for the following reasons. 

The primary and defining feature of FLAME is the accurate Trefftz ap- 
proximation; "meshlessness" is secondary. Approaches that have come to be 
known as "meshless methods" are associated with quite different approxima- 
tions. Moving least squares (MLS) techniques, the "reproducing kernel par- 
ticle method" (RKPM) [11,13,10], and, most recently, maximum-entropy 
approximations [16,17] are part and parcel of meshless techniques but are 
quite foreign to FLAME. 

Also, the similarity between least squares FLAME and the established 
meshless methods does not go very far. MLS, RKPM or maximum-entropy 
functions in meshless methods are subject to costly numerical integration 
and differentiation. In contrast, no integration or differentiation is needed 
in FLAME at aU. 

The use of the least squares fit in both MLS and FLAME is also a 
superficial similarity. In MLS, the approximation is "moving," in the sense 
that the coefficients of the relevant polynomial expansion vary from point 
to point and are found via a least squares match at the nodes. The spatial 
variation of the coefficients complicates the differentiation and integration 
that need to be carried out in the numerical procedure. In FLAME, the 
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approximation is usually non-polynomial and "static" : the coefficients are 
fixed for any given stencil. Further notes on meshlcss methods can be found 
in Section 6. 

Connections of FLAME with other classes of numerical methods such as 
variational Trefftz methods, GFEM, Discontinuous Galerkin, variational- 
difference schemes by Moskow et al. [18], discontinuous enrichment and 
others have previously been discussed in great detail [1]; see in particular 
Sections 1 & 2, Fig. 1, and references in that paper. Here I briefly comment 
on a few additional contributions, some of them recent. 

There is a well established and respectable body of work on special high- 
order schemes for various types of problems, in particular the Helmholtz 
equation [19,20,21,22,23,24]; see Harari's review [19] for further information 
and references. There is a commonality of goals but not of the methodology 
between these techniques and FLAME. 

Special difference schemes that can now be viewed as natural particu- 
lar cases of FLAME have been independently invented by various research 
groups. Mei [25] used the fundamental solutions of the Laplace equation 
to construct approximate absorbing conditions at the exterior boundary 
of the computational domain for unbounded problems. Similar ideas were 
put forward by Mittra, Boag and Leviatan [26,28,27]. Hadley [29,30] de- 
rives difference schemes for the Helmholtz equation (with applications to 
electromagnetic waveguide analysis) from the Bessel function expansion in 
free space as well as at material boundaries and corners. Very recently, 
Chiang and co-authors [31,32] developed high-order schemes for frequency 
domain analysis of 2-D photonic crystals; to do so, they carefully match first 
and higher order derivatives at curved dielectric interfaces. Also in connec- 
tion with photonic crystals, Lu et al. proposed a pseudospectral method of 
Dirichlet-to-Neumann (DtN) maps [33,34], where the field is approximated 
by cylindrical harmonics within the lattice cell and the Bloch-Floquet con- 
ditions are imposed at a set of boundary collocation points. FLAME, unlike 
DtN, relies only on local analytical approximations over a given stencil, 
which tends to be computationally more stable than approximations over 
the whole lattice cell [7]. 

The remainder of the paper is organized as follows. For completeness. 
Section 2 summarizes the established FLAME setup. Section 3 presents 
the new version of FLAME that casts the idea of Pinheiro & Webb in a 
general form. As explained in Section 4, boundary conditions in FLAME 
are easy to impose. Section 5 proves that the consistency error in FLAME 
is commensurate with the approximation accuracy of the solution by the 
FLAME bases. 

A brief review of meshless methods, in comparison with the least squares 
FLAME, is given in Section 6. A variety of numerical examples in 2D and 
3D (Section 7) include electrostatic (magnetostatic) particle interactions, 
scattering of electromagnetic (or acoustic) waves, and wave propagation in 
a photonic crystal. The examples also explore the role of the grid and stencil 
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size, of the number of approximating functions, and of the irregularity of 
the stencils. 



2 TrefFtz-FLAME Schemes: A Brief Review 

Trefftz-FLAME is a generalized finite-difference (FD) calculus that incorpo- 
rates accurate local approximations of the solution into a difference scheme. 
Conceptually, the computational domain is covered by a finite number of 
overlapping subdomains ("patches") fi'-*^, = uri^'-*, i = l,2,...?i. Each 
patch contains a stencil of a global Cartesian grid. 

Associated with each patch fi^*-* is the local approximation space 



= span{?/'W}, a = 1,2,. 



(1) 



The number m of approximating functions and the number M of stencil 
nodes may in general depend on the patch, but for simplicity of notation 
this is not explicitly indicated; i.e. we write M, m instead of M^^\ m^^\ 
Superscript (i) will occasionally be dropped in other cases as well, when it 
is clear from the context that the focus is on a particular patch (stencil) 
and there is no possibility of confusion. 

The local solution uj^''' lies in space - i.e. it is a linear combination 
of the local basis functions 



(i) 



(2) 



a=l 



In Trefftz-FLAME, the basis functions are chosen to satisfy the underlying 
differential equation, along with the interface boundary conditions. 

on stencil are 
J™. The relevant 



Clearly, the the nodal values u^^' eR of function 



linearly related to the coefficient vector c*^*^ = {cq"*} e 
transformation matrix N^'^\ 

^(0 ^ jv('\(^) 



(3) 



contains the nodal values of the basis functions on the stencil; if rk is the 
position vector of node fc, then [1,2,4] 



(4) 



A coefficient vector s*^*) G R^^ is sought to yield the difference scheme 

(5) 
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for the nodal values u^'^ of any function u]^ of form (2). Due to (3) and (5), 

For this to hold for any set of coefficients c^'' , one must have 

s(') e Null(7V(*)'^) (7) 

We shall assume that the stencil and the basis set are such that the null 
space in the definition of the scheme is of dimension one. If this happens 
not to be the case in a particular situation, the stencil and/or the basis set 
need to be changed. 

For inhomogeneous equations (i.e. with a nonzero right hand side) of the 
generic form 

Cu = /, (8) 

FLAME schemes are generated by introducing the local splitting of the 
solution: 

= + uf, (9) 

where Uq ^ is the solution for the homogeneous equation and u^'^ is a partic- 
ular solution of the inhomogeneous equation. The inhomogeneous FLAME 
scheme is [1,4] 

The coefficients of the scheme and hence the system matrix arc the same for 
the homogeneous and inhomogeneous problems. The difference is that, in 
the presence of sources in the vicinity of a given grid stencil, the right hand 
side is formed, as indicated by (10), by applying the difference operator to 
the nodal values of the particular solution Uj . 



3 Nev^f Schemes: Least Squares TrefTtz-FLAME 

This section formalizes the modification of FLAME proposed by Pinheiro 
& Webb [9] and casts it in a general form. Specific examples are presented 
in the following sections. 

One again formally considers a cover Ul^W, i = l,2,...nof the compu- 
tational domain by overlapping patches fi*-*' . Each patch again contains a 
set of nodes (stencil). We do not have to assume a regular underlying grid; 
the nodes could form fuzzy "clouds" . 

The local approximation spaces are introduced in exactly the same way 
as previously; to repeat for easy reference. 



= span{'0^')}, a = l,2,...m 



(11) 
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As before, the local solution uj^ lies in - i.e. it is a linear combination 
of the local basis functions satisfying the underlying differential equation, 
along with the interface boundary conditions: 

m 

ut' - Y.'^^^ (12) 

a=l 

At this point, the new approach parts ways with the established version 
of FLAME and involves, formally, two sets of points. The first one is n*^') 
stencil nodes iS*^*-* for which the difference scheme is being generated, and 
the other one is a reduced set 5^*^ C S^'^^ of nj*^ < n*^*^ nodes at which least 
squares matching of the solution is effected, as described below. In practice, 
the two sets will typically differ just by one extra node in 5^'-' (see below). 

For any given set of nodal values on the reduced set S{ , the linear rela- 
tionship (3) between these values and the expansion coefficients is required 
to be satisfied in the least squares (l.s.) sense rather than exactly: 

uf, ivl')c(^) (13) 

Here n\^'^ is the matrix of nodal values of the basis functions on sf^ ~ a 
matrix completely analogous to iV*^') of the previous section. The advantage 
of using the least squares match is that, as already noted, the connection 
between the set of basis functions and the stencil becomes much less rigid 
than in the established version of FLAME. 

More explicitly, the least squares relationship can be written as 

c(*) = N^^+uf^ (14) 

where N^f'^ is the generalized inverse (also known as the pseudoinverse or 
the Moore-Penrose inverse) 

iV[')+ - [n^^^ N^^y\^^^ (15) 

assuixiing that the matrix in the parentheses is nonsingular. This expression 
comes directly from the least squares methodology.^ 

The nodal values on the remaining subset of the stencil nodes Sq = 
_ g^j.g required to be exactly representable - as opposed to the 
least-squares fit - by a linear combination of the basis functions. Then the 
Euclidean vector of these nodal values is 

where iVg , directly analogous to N{ , contains the values of the basis 

{i) 

functions on the nodes of S'q . For a homogeneous differential equation 

^ It is well known that matrix inversion is hardly ever needed in practical 
computation. In particular, the pseudoinverse can be computed using the QR 
method rather than by explicit matrix inversion. 



8 



Tsukerman 



(zero right hand side), we are again looking for a coefficient vector s*-*^ of 
the difference scheme 

s(')V^) = (17) 

or, equivalently due to (16) 



Inl 

where /„i is the rii x ni identity matrix. A natural partitioning of the 
difference scheme s*^*'-^ — (sq''"^, s^i^"^) is assumed, with the coefficients 
Sq^'^ , Si^"^ corresponding to the nodes in Sq^ and s[^\ respectively. 

Since the vector of nodal values u^g\ is unknown and in a sense quasi- 
arbitrary, we require that the above relationship be satisfied for any such 
vector. This immediately gives 



s 



W G Null(^o'^^i*'^^) ^ Null(ArW+^ArW^, (19) 



This is our main expression generalizing the null space condition (7) of the 
established version of Trefftz-FLAME. We shall always assume that the 
stencil and the basis arc set up in such a way that the null space in the 
definition of the scheme is of dimension one. Otherwise the basis and/or the 
stencil need to be modified. 

It immediately follows from (19) that the scheme can be explicitly writ- 
ten as 

= (4'^"^, -.«^iV«7V«+) (20) 

(i) . . 

Typically Sq will contain only one node, in which case, since the null space 
is defined up to an arbitrary factor, one can set Sq'^ = 1 and obtain 

^ (1, -iV«iv(*)+) (21) 

This is equivalent to Pinheiro & Webb's algorithm that they used for elec- 
tromagnetic wave scattering in 3D [9]. 

The rectangular matrix on the right of (19) has the number of columns 
equal to the stencil size ti*-*-* and the number of rows equal to the number 
of nodes n^*^ in the substencil Si{i). Ordinarily, such a matrix will have 
a one-dimensional nuUspace (leading to a unique FLAME scheme) if the 
numbers of columns and rows differ by one - that is, if n^*^ = n*-*^ — 1 and 
S'o(i) has exactly one node, as stipulated above. 

For convenience. Figs. 1-4 illustrate the dimensions and composition 
of the matrices involved in the key equations from (14) to (19). In this 
illustration, there are four approximating functions and seven stencil nodes, 
five of which are the "least squares" matching nodes of ^i. Although one 
node rather than two would be typical for Sq, two nodes are assumed, to 
demonstrate the general principle. 
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Fig. 1 Matrix dimensions in (14) for tIj*' = 5, m : 



^51 



(0 



to 



Fig. 2 Matrix dimensions in (16) for ri{^ = 5, Ji'*' = 7, m = 4. 
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Fig. 3 Matrix dimensions and structure in (18) for nj^'' = 5, n'*' = 7. 
4 Boundary Conditions in FLAME 

Boundary conditions at the material boundaries are built into the TrefFtz- 
FLAME bases and therefore need no other treatment, in sharp contrast with 
the methods where smooth functions not obeying the boundary conditions 
are employed (e.g. MLS or radial basis functions). Exterior boundary condi- 
tions in FLAME do not pose any problem either. Since all approximations 
are local, approximations near exterior boundaries are completely decou- 
pled from approximations elsewhere. In fact, any traditional non-FLAME 
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Fig. 4 Matrix dimensions and structure in (19) for Vi{^ = 5, n'*' — 7 



schemes, including perfectly matched layers (PML), can be used at the 
boundaries. One could also opt for FLAME schemes at the boundary - the 
FLAME PML [2,4] being one example. 

Natural boundary conditions (Neumann and the like) can either be han- 
dled in the established ways by traditional schemes or, alternatively, using 
FLAME with a subset of Trefftz functions satisfying the natural conditions. 

In the customary algorithmic implementation of Dirichlet conditions, 
one initially disregards the distinction between the nodes at the exterior 
boundary and the inner nodes. Then, after the system matrix has been 
assembled, the usual adjustments are made: terms of the form ayitj, where 
Uj is a given Dirichlet value and is a matrix entry, are moved over to the 
right hand side. The entry in the right hand side is modified accordingly 
and the matrix entry 0,^ is then set to zero. 

Even "fuzzy" exterior boundaries, where the diffuse set of nodes does 
not fall exactly on a predefined surface, could be handled in the same fash- 
ion. There may be some redundancy if the fuzzy boundary layer is slightly 
"thicker" than necessary in some places. It is trivial but optional to elim- 
inate the redundant nodes (for example, the ones connected only to other 
Dirichlet nodes) from the system. 



5 Consistency Error of FLAME 

Analysis of the consistency error of least squares FLAME proceeds along 
the same lines as in the original version of FLAME [1,2], with some mod- 
ifications. Here we shall only consider equations with the zero right hand 
side in a given patch (i.e. over a given grid stencil); the treatment of equa- 
tions with nonzero r.h.s. is exactly the same as in the established versions 
of FLAME [1,2]. 

The consistency error of scheme (17), (21) is, by definition, obtained by 
substituting the nodal values of the exact solution u* into (17). Let ta{K) 
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be the approximation error of the exact solution u* within a patch rj^'': 



ea{h) = min - ^ c^^W (22) 
Equivalently, there exists a coefficient vector c*-*^ £ M™ such that 

U* - \Moo -^a{h) (23) 

^ — 'a— 1 

within the patch. For the nodal values, one then has due to (16) 

AA(*)u* = Af(*)c(*) + Ti (24) 

where Af^^^u* and ry = Jsf^^^r] are the Euclidean vectors of nodal values of 
the exact solution and of 77 over stencil i, respectively. A^''-* is (as always) 
the matrix of nodal values of the basis functions. Due to (23), 



ll^^lloo < ea(/i) 



In particular, for the nodal values on the "least squares part" S^'^ of the 
stencil, we have 

^*(*) ^ ^ (25) 

Left-multiplying this equation with N^^"'" and then expressing c^'^ yields 

c« = (iV«^iV«)"'7V«^(^,;« _^^) EE iV«+(^.;« ~^^) (26) 

assuming that the matrix in the parentheses is nonsingular. The ith com- 
ponent of the consistency error for scheme (21) is, by deffirition,'^ 



\f-ci{h)\ 



= (7V^^^c(') + 770) - TVf ivf ^+(ivf ^c(') + 7^^) 
Now substituting c'*^ from (26), we have 

-iv(%«+^J = |r7o-iV(^)iv(^)+^J < |77o| + ||A^^*^iVp|| llr/JI (27) 

where the term in the square brackets has in the end vanished because 
= for any pseudoinverse matrix . As a result of the general- 
ized inverse falling out of the accuracy estimate, it can be seen - somewhat 
counterintuitively - that the least squares fit does not reduce the accu- 
racy as compared to the exact nodal match. The consistency error is still 
governed primarily by the approximation error. 



Analysis for the slightly more general scheme (20) is similar. 
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6 Trefftz-FLAME vs. Meshless Methods 

Meshless methods have been developed and studied very extensively over 
the last two decades, and excellent contributions and reviews [11,12,13,14, 
15], including a detailed recent one [10], are available. This section is limited 
to a brief summary of such methods, with the emphasis on some similarities 
with FLAME but - more importantly - on the differences. 

Meshless techniques have two key advantages. One is, by definition, the 
absence of complex meshes that may be difficult to generate, especially 
in 3D, - although regular auxiliary grids are still needed to compute the 
numerical quadratures in weighted residual and Galerkin methods. 

The second advantage is high-order approximation, at least in the re- 
gions where the solution is smooth. This can be accomplished by the Re- 
producing kernel particle method (RKPM), Moving least squares (MLS) 
or similar approximations, as well as by the recently developed maximum- 
entropy functions [16,17]. 

However, meshless methods have problems that are, in a certain sense, 
an extension of their advantages. Approximation by smooth functions is 
accurate in homogeneous regions but breaks down at interface boundaries 
where the solution or its derivatives have jumps. The MLS, RKPM and simi- 
lar approximating functions and especially their derivatives are too complex 
to be differentiated or integrated analytically. Therefore one has to resort to 
numerical integration and differentiation, which is in general quite costly. 
Not only interface boundary conditions but also the Dirichlet conditions 
on the exterior boundary of the domain are difficult to impose numerically 
because the basis function do not satisfy the Kronecker delta property at 
the nodes. The maximum-entropy functions of [16,17] satisfy a "weak Kro- 
necker delta" property due to their rapid decay; however, to compute these 
functions, one has to solve a nonlinear optimization problem. 

These impediments and ways of getting around them have been the 
subject of much research that is still ongoing [10]. The difficulties can be 
alleviated but not removed. For example, if collocation instead of weighted 
residual / Galerkin methods is used, numerical integration becomes unneces- 
sary but, as an unpleasant trade-off, higher-order derivatives, commensurate 
with the order of the differential operator, of the basis functions have to be 
calculated. (In the weighted residual methods, half of the derivatives are 
usually moved over from the basis functions to the test functions, thereby 
reducing the required order of differentiation. For collocation, this avenue 
is not available.) 

Interface boundary conditions in meshless methods can be treated using 
Lagrange multipliers or by introducing additional "enrichment" functions. 
(The latter can also be used to represent boundary layers, singularities, etc.) 
This, however, increases the computational and algorithmic complexity. La- 
grange multipliers not only constitute additional unknowns in the system 
but may also adversely affect the algebraic properties of the system and may 
lead to ill-conditioning. The enrichment functions may have to be mollified 
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by partition of unity factors (to make their support small) and then are also 
subject to numerical integration and differentiation. 

In contrast, TrefFtz-FLAME basis functions satisfy the underlying equa- 
tion and the boundary conditions by construction and therefore do not 
require any additional constraints to be imposed. In essence, all FLAME 
functions are "enrichment" functions in the sense outlined above; all of them 
are physical and reflect the local behavior of the solution. No degrees of free- 
dom are wasted on approximating broader classes of functions not directly 
relevant to the problem under consideration. 

In short, the difflculties that traditional meshless methods have in deal- 
ing with interface and essential conditions are real, but ways have been 
developed to overcome that. FLAME does not have these difficulties to 
begin with, so there is nothing to overcome. The price to pay for this ad- 
vantage is the need to pre-define accurate local analytical or semi-analytical 
bases. The complexity of this task depends on the class of problems, making 
FLAME quite suitable for some classes and less so for the others. Exam- 
ples of problems where the new version of FLAME works well are presented 
in Section 7; see also [1,2,4,5] for many other examples of the established 
version of FLAME. 

7 Computational Examples of TrefFtz-FLAME on Irregular Grids 

7.1 2D Electromagnetic Wave Scattering 

Let us start by applying least squares FLAME to the classical problem of 
electromagnetic scattering from a dielectric cylinder. Non- ideal conductors 
can be treated in the same manner as dielectrics if complex permittivity is 
used. The wave is assumed monochromatic 

In the 2D case, where the material parameters and fields are independent 
of one coordinate (say, z), electromagnetic waves can be decomposed into 
two modes. In the i?-mode (known as the s-mode in optics and frequently, 
but not always, also referred to as the TM mode) the electric field has 
only one component E ~ Ez, whereas the magnetic field H has x and y 
components, with = 0. Similarly, for the iJ-mode (also known as the p- 
or TE mode) H ^ H^, E^ = 0. These modes satisfy the following equations 
that are easy to derive from Maxwell's system: 

V -fi-^VE + uj^eE = (28) 

V • e-^\7H + uj^fiH = (29) 

For numerical testing, let us consider a dielectric cylinder of radius rcyi and 
dielectric permittivity Ccyi- The host medium is assumed to have the relative 
dielectric constant of one, eout = 1- The incident field is a plane wave with 
a wave vector k, 

Einc = Eoexp(ik-r) 
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under the exp{—iujt) complex phasor convention. To complete the formu- 
lation, the standard Sommerfeld radiation conditions can be imposed on 
the scattered field Esc = E — Einc- Numerically, well established techniques 
such as PML or absorbing boundary conditions can be applied on the exte- 
rior boundary. Since our focus is on the construction of accurate difference 
schemes in the interior, absorbing conditions are a tangential issue, and in 
the numerical tests the exact solution was imposed as a Dirichlet condition. 
(A textbook solution for a cylindrical scatterer can be found e.g. in [35].) 
When Dirichlet conditions are imposed, it is tacitly assumed that w is not 
an eigenvalue of the corresponding operators V • /i~^V or V • e~^V. 

The computational domain is normalized to the unit square [—0.5,0.5] x 
[—0.5,0.5], In all numerical experiments, the scatterer is centered at the 
origin and has the relative permittivity of Ccyi = 9 and radius Tcyi = 0.25. 

Irregular grids are generated by randomly displacing the nodes of a 
regular Cartesian grid: 

Xi ^ Xi + JrjxiK (30) 

where Xi is the coordinate of the ith node of the regular grid, is the 
grid size in the x-direction, rjxi is a random number uniformly distributed 
within (—0.5,0.5). Parameter / controls the "fuzziness" of the distribution 
and is chosen between and 0.5, to maintain the majority of the topological 
connections between each node and its nearest neighbors. Displacements in 
the y direction are generated in a completely similar way. Fig. 5 gives an 
example of a 12 x 12 irregular grid. 
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Construction of the TrcfFtz bases has been elaborated upon in the pre- 
vious publications [1,2,4,5]. Briefly, the basis functions are chosen as cylin- 
drical harmonics: 



where J„ is the Bessel function, Hn is the Hankel function of the first 
kind"*, and a„, b„ are coefficients found via the standard conditions on the 
boundary of the cylinder. The FLAME basis contains a finite number of 
harmonics from the monopole (n — 0) up to some highest order rimax- 

The relative numerical errors for the i^-mode computed on a 30 x 30 
FLAME grid, with a varying number of basis functions, stencil nodes and 
fuzziness, are reported in Table 1. The stencil of size m for any given node 
is formally defined as the set of m nodes closest to that node (including the 
node itself). The relative errors were computed as 



The overall high accuracy of FLAME is evident from the table. Five to ten 
correct digits in the value of the field are routinely obtained, even though 
the grid is irregular, relatively coarse and does not represent the geometry 
of the problem in any way. Rather, the geometric and physical information 
is built into the FLAME bases. 

While the high accuracy was to be expected, there are two surprises in 
the results of Table 1. First, the fuzziness of the grid does not significantly 
affect the accuracy - if anything, the errors on fuzzier grids (/ = 0.2, / = 
0.4) tend to be a little lower than on the regular one (/ = 0). No "physical" 
reason for that is apparent, especially in light of the fact that the situation 
in 3D is different (see the 3D example below). 

The second surprise is that, while the accuracy improves as expected 
when the number of multipoles rimax is increased from three to four to five, 
the relative error goes up from ~ 10~^^ to ~ 10~^° when the 6th order 
multipoles are included. One possible circumstance that might contribute 
to this was noted by Cajko [36]: the coefficients of high-order schemes must 
themselves be computed with very high precision, to avoid an asymptotic 
loss of accuracy for fine grids. 

7.2 Wave Propagation in Photonic Crystals 

Photonic crystals (PhC) are artificial periodic structures that exhibit very 
peculiar characteristics of wave propagation (e.g. the photonic bandgap [37, 




^^'^ = a„J„(fccyir) exp(in0), r < ro 
[6„ J„(A:air?') + H^^^ (fcairf )] exp{in(f>), r > To 



ll*FLAME ~ licxactll 



(31) 



1 1 Inexact 1 1 



* Hankel functions of the second kind should be used if the exp(itjt) phasor 
convention common in electrical engineering is adopted. 
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7^max (highest harmonic 
in the FLAME basis) 


Stencil size 


"Fuzziness factor" / 


Numerical error 


3 


9 





l.lE-06 


4 


13 





2.3E-09 


5 


21 





1.29E-11 


6 


21 





4.5E-10 


3 


9 


0.2 


l.lE-06 


4 


13 


0.2 


6.2E-11 


5 


21 


0.2 


1.5E-11 


6 


21 


0.2 


3.6E-10 


3 


9 


0.4 


8.7E-07 


4 


13 


0.4 


7.3E-11 


5 


21 


0.4 


1.3E-11 


6 


21 


0.4 


4.7E-10 



Table 1 Relative numerical errors in the nodal values of the electric field for 
least squares FLAME schemes with varying bases and stencils on a 30 x 30 grid. 
Scattering from a dielectric cylinder with tcyi = 9, rcyi = 0.25. 



38,39,40]) and are finding various applications in lightwave technology [41. 
42,2]. Our goal here is limited to using PhC as a demonstration example 
for the new version of FLAME. 

For convenience of comparison, consider the same numerical example as 
in [2,4]: a photonic crystal due to Fujisawa & Koshiba [43], with a square 
lattice of cylindrical coaxial dielectric rods and a bended waveguide obtained 
by eliminating a few of these rods (Fig. 6). The E- and iJ- modes of the wave 
are governed by equations (28) and (29), respectively. 

As in the previous tests [2,4], the dielectric constant is set to erod = 
9 for the rods (index of refraction n = 3) and to unity for the outside 
medium. The radius of the cylinders and the wavenumber are normalized 
to unity; the air gap between the neighboring rods is equal to their radius. 
The field distribution is shown in Fig. 6 for illustration. Simplified boundary 
conditions are set equal to an externally applied plane wave. This is adequate 
for the demonstration example; conditions at the ports of the guide are a 
separate issue, and a way to handle them in FLAME was developed by 
Pinheiro & Webb [44]. 

An excellent agreement between the established version of Trefftz-FL AME 
and independent FEM simulations has been previously demonstrated [2,4]; 
in fact, numerical evidence was convincing that the numerical error was 
in this case orders of magnitude lower in FLAME than in FEM with a 
comparable number of unknowns. 

Our goal here is to compare the least squares FLAME on fuzzy node 
clouds with the previous Trefftz-FL AME on regular Cartesian grids, thereby 
exploring the influence of the "fuzziness" on the numerical accuracy. The 
"fuzzy" node coordinates Xi, jji are generated as in (30). Fig. 7 shows an 
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Fig. 6 The imaginary part of the electric field in the photonic crystal waveguide 
bend. The real part looks qualitatively similar. (Reprinted by permission from [4] 
©2005 IEEE.) 



example of an irregular 60 x 60 FLAME grid with the "fuzzincss factor" 
/ = 0.4. 

Tliis grid was used to produce the least squares FLAME results in Fig. 8, 
where the distribution of the electric field along the midline of the crystal 
(y = 0) is shown for the £'-mode. The result obtained with the new version 
of FLAME coincides with the various results obtained previously with FEM 
and FLAME, despite the quite irregular distribution of nodes in the new 
test (see Fig. 7). 

Note that for the 60 x 60 grid there are about 12.5 points per wavelength 
(ppw) in the air but only 4.2 ppw in the rods, and yet the FLAME results are 
very accurate even on an irregular grid, due to the Trefftz approximation. 
Alternative methods that rely on generic polynomial approximations would 
require a substantially higher number of ppw to achieve the same accuracy. 

The FLAME field values plotted in Fig. 8 were produced by interpo- 
lation, similar to the way described in [2], pp. 290-291, for the previous 
version of FLAME. Namely, for a given set of nodal values on a stencil, one 
finds the Trefftz expansion coefficients c*^*^ by solving the least squares prob- 
lem either on the substencil Si (13) or, alternatively, on the whole stencil. 
The interpolation is then effected by the Trefftz expansion (2). 
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60 X 60 fuzzy grid with / = 0.4 




X 



Fig. 7 A 60 X 60 irregular FLAME grid in the photonic crystal example. "Fuzzi- 
ness factor" / = 0.4. 



Fujisawa & Koshiba crystal 
Real part of the field along midline (y = 0) 
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• real(E) FLAME 150x150 

* real(E) FLAME 50x50 
° real(E) FLAME 30x30 



— real(E)FEMLAB 154461 
dofp=2 

— real(E) FEMLAB adaptive 
1 07604 dofp=2 

^real(E) FEMLAB 33679 dof 
p=2 

□ 17-ptl.s.-FLAME, 60x60, 11 
b.f. 



Fig. 8 Field distribution in the Fujisawa-Koshiba photonic crystal along the 
central line y = 0. Solution by the new version of FLAME coincides with the 
previous FEM and FLAME results, despite the quite irregular distribution of 
nodes in the new test. For testing, least-squares FLAME was applied on 17-point 
stencils with 11 basis functions. 
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7.5 Electrostatic Particle Interactions in 3D 

7.3.1 Formulation and setup. The 3D numerical example, chosen to be 
simple enough for illustration but extendable to more complex cases, in- 
volves two spherical dielectric particles in a uniform external electrostatic 
field. 

A uniform field induces only the dipole mode in a single particle [35] , but 
for two or more particles higher-order multipoles do arise. A semi-analytical 
solution can be obtained by the well known multipole-multicenter method 
(mmc) [45,46,2], where the potential is sought as a superposition of multi- 
pole expansions centered on each of the particles. To impose the boundary 
conditions on any given particle, one needs to translate the multipole ex- 
pansions from all other particles to that particle; this is accomplished by 
the standard multipole translation formulas (e.g. [45]). Since the details of 
this procedure are easily found in the literature and are tangential to the 
subject matter of this paper, they are omitted here. 

The electrostatic equation in terms of the total electric potential u is^ 

V ■ eVit = 0; u(r) — Ucxt(r) — > as r ^ cx) (32) 

where Ucxt(r) is the applied (external) potential, which for a uniform field is 
a linear function of coordinates. If understood in the sense of distributions, 
(32) implicitly includes the standard conditions on interface boundaries: the 
continuity of the potential and the normal component of the electric flux 
density D — eE. The div-grad equation (32) can also be rewritten in terms 
of the "scattered" potential Ug = u — Ucxt, but there is no particular need 
to do so here. 

The construction of FLAME bases is analogous to the cylindrical case 
and involves spherical harmonics defined inside and outside a given particle 
and matched at its interface. Details can be found in [4,2]. It is possible to 
increase the accuracy of FLAME by using more complex mmc bases over 
several neighboring particles as done in [8]; however, in this paper single- 
center basis functions are used for the sake of practical simplicity. 

The mmc expansion with a large number of harmonics (64 terms per 
particle in the expansion) is used to compute a quasi-analytical solution for 
reference and error analysis. 

The two spherical particles in our 3D example have radii ri = 1 and r2 = 
1.25 and are centered on the z axis at points zi — —2.5, Z2 ~ 2, respectively. 
The lack of geometric symmetry is deliberate, to avoid any computationally 
favorable artifacts due to symmetry. The dielectric constants inside and 
outside the particles in all numerical experiments were chosen as e;,! = 5, 
£out = Ij respectively. A cubic computational domain of size idom = 10 
centered at the origin was used. 

^ u is used instead of the somewhat more conventional <^ to avoid confusion with 
the polar coordinate. 
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l.OOE-06 1 

0.1 grid size 1 

Fig. 9 Relative error in potential, regular grids. 19-point least-squares FLAME 
with 9 basis functions. Test for two spherical particles. To evaluate the numerical 
error, the mmc solution with 64 terms per particle was taken as quasi-exact. The 
dashed line is a visual aid indicating the O(h^) convergence rate. 



7.3.2 Regular Cartesian grids. The first set of tests is performed on regular 
Cartesian grids. The relative simplicity of this example notwithstanding, 
note that the regular grids do not carry any accurate geometric information 
about the shapes of the particles; that information is implicitly imbedded 
in the FLAME approximating functions. In contrast with the established 
version of FLAME, now the stencil size is to a large degree independent of 
the number of basis functions. 

Fig. 9 illustrates convergence of 19-point least squares FLAME with nine 
basis functions. The relative error plotted in the figure was calculated as 

r _ ll:?iFLAME ^ limmcll fnn\ 



where Mflame ^'^d l^mmc Euclidean vectors of the nodal values of FLAME 
and mmc potentials, respectively; the norms are EucUdean. 

The dashed line in Fig. 9 is a visual aid indicating the 0{h^) slope and 
demonstrating quadratic convergence, consistent with the number of basis 
functions (nine) used. 

7.5.5 Irregular stencils. A fragment of an irregular distribution of nodes 
near one of the two particles is shown in Fig. 10. The random displace- 
ments of the nodes are generated in the same manner as in the previous 
examples, except that now these displacements are also applied in the z 
direction in addition to x and y. The exact potential was imposed as the 
Dirichlet boundary condition for testing, to avoid any extraneous errors due 
to domain truncation. 

Table 2 demonstrates the effect of "fuzziness" on the numerical accu- 
racy. For all combinations of basis functions and stencils in the table, the 
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Fig. 10 An irregular set of nodes for FLAME in the vicinity of a particle. The 
nodes are randomly displaced with respect to the regular Cartesian grid 40 x 40 x 
40, with the "fuzziness" (see text) of f — 0.4. Nodes inside the particle not visible 
but are distributed in a similar manner and with the same density. 



Number of 
basis functions 


Stencil size 


/ = o 


/ = 0.1 


1 = 0-2 


/ = 0.4 


6 


9 


7.9E-04 


8.05E-05 


1.22E-04 


2.69E-04 


6 


19 


4.72E-06 


5.31E-06 


1.06E-05 


2.49E-05 


6 


27 


9.69E-06 


9.17E-06 


1.67E-05 


1.74E-05 


9 


19 


2.13E-06 


2.00E-06 


3.10E-06 


9.17E-06 


9 


27 


1.14E-05 


1.12E-05 


9.84E-06 


6.68E-06 



Table 2 Relative errors in the nodal potential for FLAME on irregular grids. The 
reference quasi-analytical solution obtained by the multipole-multicenter (mmc) 
expansion with 64 terms per particle. The 3D test problem with two dielectric 
particles with the dielectric permittivity 5. 41^ nodes in FLAME. 



numerical accuracy deteriorates somewhat as the "fuzziness" / increases. 
This is different from the previous 2D example where the fuzziness was not 
a significant factor. 

The influence of the stencil size for a fixed number of basis functions is 
less straightforward. For six basis functions, the accuracy tends to improve 
when the stencil size changes from 9 to 19 but deteriorates somewhat as the 
stencil expands further. 

A plausible qualitative explanation is as follows. As the number of stencil 
nodes increases, so does the robustness of the analytical approximation in 
FLAME, because more information (in the form of the nodal values) is being 
utilized. On the other hand, as the geometric size diam(rj('^) of the stencil 
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"patch" increases, the approximation accuracy, proportional to diani(r2'^*')'', 
deteriorates {p is the order of approximation). These two conflicting trends 
lead to the nontrivial dependence of the accuracy on the stencil size. It 
would be worthwhile to analyze this rigorously in the future. 

7.3.4 Random nodes. Further pushing the envelope, I have tested a com- 
pletely random distribution of nodes (as opposed to random displacements 
around a regular lattice), see Fig. 11. Clearly, this can be viewed as the 
worst case scenario, and much better results can be expected if sensible 
adaptive procedures [8] are adopted in the future. 

The initial test is performed with one particle only. In this case, since 
the exact solution contains only the dipole harmonic and is thus included 
in the FLAME basis, one can expect that FLAME will produce the exact 
solution up to the round-off error. This has been confirmed experimentally. 
For example, for 5000 FLAME nodes and the 27-point FLAME scheme with 
16 basis functions, the relative error in the potential at the nodes was on 
the order of 10""'^'^, despite the fact that the distribution of the nodes is 
random and hence carries no geometric information whatsoever. 

Note that 5000 nodes correspond to only ~ 17 points per coordinate, 
with the average separation distance of ft, idom/17 ~ 0.59, which is 
comparable with the radius of the particle. Clearly, classical FD schemes 
would require much finer meshes for high accuracy (and still would not be 
able to deliver machine precision). 

Let us now return to the test problem with two particles, as described 
above. Fig. 11 presents the general setup of the problem solved with a 
random distribution of 8000 nodes (only 1600 shown in the picture); the 
simulation results are reported below. Fig. 12 shows projections of the same 
nodes and of the spheres on the xy-plane. 

Table 3 shows the dependence of the numerical error on the number of 
random nodes. The errors are still quite small despite the highly suboptimal 
distribution of nodes. At the same time, as could be expected, the errors are 
higher than for regular grids and, moreover, convergence of the solution with 
respect to the number of nodes is not monotone. This is almost certainly 
attributable to "poor quality" stencils, although precise a priori measures 
of quality in FLAME are still to be developed; see further comments on 
that below. 

One may envision that in future applications of the new version of 
FLAME the nodal distributions will be optimized and generated adaptively, 
based on appropriate a posteriori error estimates. Such distributions can 
then be expected to produce more accurate solutions than random, semi- 
random or even regular sets of modes would. 

A natural a posteriori error indicator'' in FLAME is the discrepancy 
between the values of the solution over two overlapping patches; this indi- 

® As in FEM, one may want to distinguish estimators that produce some quan- 
titative measure of the error from qualitative indicators of regions where grid 
refinement is desirable. 
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Fig. 11 Two dielectric spheres in an external field. An example of least squares 
FLAME with 8,000 randomly generated nodes (only 1,600 shown for visual clar- 
ity). Nodes inside the particle not visible but are distributed with the same density. 



Number of nodes 


Effective "grid size" 


Relative error 


10000 


0.464 


1.94E-04 


20000 


0.368 


3.66E-04 


30000 


0.322 


2.76E-05 


40000 


0.292 


3.16E-05 


60000 


0.255 


2.11E-05 



Table 3 Relative errors in the nodal potential for FLAME on fully random sets 
of nodes. The reference quasi-analytical solution is obtained by the multipole- 
multicenter (mmc) expansion with 64 terms per particle. 3D test problem with 
two dielectric particles with the dielectric permittivity 5. 19-point stencils with 9 
basis functions. 



cator has already been shown to produce quite reasonable results, leading 
to grid refinement around the gaps between two neighboring particles [8, 
47]. It should certainly be possible to develop other, and better, indicators 
in the future. 

Likewise - and also in parallel with FEM where similar issues have been 
studied for several decades ~ good a priori error estimates indicating the 
"quality" of irregular stencils are needed. In FEM, there are several geo- 
metric conditions that characterize the quality of elements ( "quality" being 
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Fig. 12 Same as in the previous figure, but projections of tlie nodes and of tlie 
two particles on tfie xy-plane sliown. Note tliat tlie exterior Dirichlet boundary 
in this example is "fuzzy" . 

ultimately understood as approximation accuracy): for example, the mini- 
mum or maximum angle conditions, the radius of the inscribed sphere vs. 
the element diameter, etc. 

In comparison with FEM, irregular stencils (especially large ones) are 
difficult to characterize in pure geometric terms. Therefore particularly rel- 
evant are the algebraic conditions. In FEM, such conditions are related to 
the properties of the affine transformation of a given element to the canon- 
ical "master" element [48,49] or. in less traditional analysis (see [2], §3.14), 
to the maximum eigenvalue of the element stiffness matrix or the minimum 
singular value of the edge shape matrix [52,53,51,50]. One conjecture that 
may form an interesting subject for future research is that for least squares 
FLAME the minimum singular value of the nodal matrices such as TV (4) 
will play a similar role. 



8 Conclusion 

FLAME replaces standard Taylor expansions of classical finite difference 
methods with local approximations of the solution by functions satisfying 
the underlying differential equation and interface boundary conditions. As 
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a result, it is not unusual for the consistency error of FLAME to be orders 
of magnitude lower than that of conventional schemes, due to the high 
accuracy of the Trefftz approximations employed in FLAME. 

Previously one constraint in FLAME has been that the number of ap- 
proximating functions had to be closely linked to the number of nodes in 
the grid stencil. The present paper, by expanding on the idea of Pinheiro 
& Webb [9], removes this constraint: while Trefftz functions still provide 
very accurate approximation, the numbers of functions and nodes are now 
decoupled. This is accomplished by using least-squares matching of the ap- 
proximate solution at the stencil nodes. Consequently, FLAME schemes can 
now be generated on irregular stencils with the number of nodes substan- 
tially greater than the number of approximating functions. This extends the 
practical applicability of the method. 

Counterintuitively, the use of the least squares fit instead of the exact 
nodal match does not reduce the numerical accuracy; on the contrary, the 
method tends to become more robust. The consistency error is still governed 
primarily by the approximation accuracy of the solution with the Trefftz- 
FLAME basis. 

In contrast with meshless methods and other similar techniques, the 
treatment of interface and exterior boundary conditions in FLAME is simple 
and natural. Trefftz-FLAME basis functions satisfy the underlying equation 
and the boundary conditions by construction and therefore do not require 
any additional constraints to be imposed. 

In summary, the new version of FLAME is a promising technique that 
combines the high accuracy of Trefftz approximations with robust operation 
on irregular stencils. Once a set of basis functions is established, FLAME 
is exceptionally easy to implement in computer codes, as it does not re- 
quire mesh generation, numerical integration or differentiation. The ability 
of the method to handle irregular stencils, as well as the a posteriori error 
indicators inherent in FLAME, bodes well for the development of adaptive 
algorithms in the future. 

The price to pay for these advantages is the need to derive the local 
Trefftz sets. Expanding the library of canonical solutions (e.g. elliptic prob- 
lems / wave scattering; (piecewise) planar / cylindrical / spherical bound- 
aries, etc.) may lead to the commensurate expansion of the applicability of 
FLAME. On the other hand, problems where local solutions are unavailable 
or difficult to obtain may remain out of reach for FLAME until additional 
ideas are put forward. 

An interesting subject for future research is a priori and a posteriori 
error estimates in FLAME. Some considerations are already presented in 
this paper. 

In comparison with the previous version of Trefftz-FLAME, least squares 
FLAME is particularly promising in cases where an intrinsic disparity be- 
tween the natural sizes of the bases and stencils exists. Examples include 
finite difference time domain simulation of wave propagation, as well as 
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electromagnetic vector problems where this version of FLAME has already 
been explored [9]. 
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